##########################################################################################
#Replication code for empirical analysis in Becher and Christiansen,  
#"Dissolution Threats and Legislative Bargaining", forthcoming at AJPS

#Written by Michael Becher (michael.becher@uni-konstanz.de)

#Version: May 27, 2014
##########################################################################################

#Required data file: dissthreats_data.Rda
#Required packages: Zelig, version 3.5.5
#The simulation of the quantities of interest reported in Table 2 and Figure 3 of the paper were produced in R using the Zelig package, version 3.5.5 (Kosuke Imai, Gary King, and Olivia Lau, "Zelig: Everyone's Statistical Software," http://gking.harvard.edu/zelig). Using the exact code in this file requires installing Zelig version 3.5.5, and de-installing newer versions, as Zelig keeps changing.
#Download version 3.5.5 from CRAN (http://cran.r-project.org/src/contrib/Archive/Zelig/Zelig_3.5.5.tar.gz) and save on computer
#To remove newer version of Zelig, type
remove.packages("Zelig",)
#Install previous version from source - specify relevant folder(s) if not saved in R directory
install.packages('Zelig_3.5.5.tar.gz', repos = NULL, type="source")
###


#Load package
library(Zelig)

#Set working directory that contains file with replication data using setwd()

#Read data set
Data <- readRDS("dissthreats_data.Rda")


##########################################################
#Figure 2: Dissolution Threats in Denmark
############################################################
#Aggregate dissolution threats by year
threats.total.year <- tapply(Data$nthreat, Data$year, FUN=sum) 
#Descriptive barplot
barplot(threats.total.year, 1974:2011,
ylab = "Number of Dissolution Threats",
xlab = "",
ylim = c(0,4),
las = 2, #x-axix labels perpendicular to axis
#label bars - for readability, print only label for every other year
names.arg = c("1974", "" ,"1976","" ,"1978","" ,"1980", "","1982","" ,"1984", "","1986","" , "1988", "","1990","" ,"1992", "","1994", "","1996", "","1998","" ,"2000","" ,"2002", "","2004","" ,"2006","" ,"2008", "","2010", ""), 
cex.names = 1
)


#############################################################################
#Table 1: The Conditional Effect of Public Opinion on Dissolution Threats: Logit Estimates
#############################################################################


m1 <- zelig(threat ~ po + cabseats + expir + dpmcp , model = "logit", data = Data)
m2 <- zelig(threat ~ po*cabseats*expir + dpmcp, model = "logit", data = Data)
m3 <- zelig(threat ~ po*cabseats*expir + dpmcp + pmparty, model = "logit", data = Den)
m4 <- zelig(threat ~ po*cabseats*expir + dpmcp + year, model = "logit", data = Data)
m5 <- zelig(threat ~ po*cabseats*expir + dpmcp + threat_lag, model = "logit", data = Data)
m6 <- zelig(threat ~ po*cabseats*expir + dpmcp + pmparty + year + threat_lag, model = "logit", data = Data)

#Print coefficient estimates
summary(m1)
summary(m2)
summary(m3)
summary(m4)
summary(m5)
summary(m6)



#############################################################################
#Table 2: Predicted Probability of Dissolution Threat by Prime Minister
#############################################################################

#We use quasi-Bayesian simulation to caluculate the predicted probability of dissolution threat given low/high seat share and low/high popularity of the government

##First, we obtain coefficient estimates (repeated from above) using model 2 from Table 1
m2 <- zelig(threat ~ po*cabseats*expir + dpmcp, model = "logit", data = Data)
##Second, we set the variables to the counterfactual values of interest
#Public opinion: low (25th percentile) or high (75th percentile)
po_low <- 0.33
po_hi <- 0.40
#Cabinet seats: low (25th percentile) or high (75th percentile)
cab_low <- 0.33
cab_hi <- 0.43
#Thus, there are four counterfactual scenarios represented by four distinct vectors:
x.2a <- setx(m2, po = po_hi, cabseats = cab_low, expir = 14)
x.2b <- setx(m2, po = po_low, cabseats = cab_low,expir = 14)
x.2c <- setx(m2, po = po_hi, cabseats = cab_hi, expir = 14)
x.2d <- setx(m2, po = po_low, cabseats = cab_hi, expir = 14)
#Note that expir is always set to 14, which corresponds to third quartile

##Third, simulation. Draw coefficient estimates from multivariate normal distribution (where mean corresponds to those estimated above) and caluclate predicted value for each scenario and their differences. Do this 10,000 times and store the results. Then calculate the mean and 95 percent confidence intervals. The results are stored in a matrix called Tab.PP
Tab.Mean <- matrix(ncol = 3, nrow =3, dimnames = list(c("Weak gov", "Strong gov", "Difference"), c("Popular gov", "Unpopular gov", "Difference") ) )
Tab.CI <- matrix(ncol = 6, nrow =3, dimnames = list(c("Weak gov", "Strong gov", "Difference"), c("Popular gov LB", "Popular gov UB",  "Unpopular gov LB", "Unpopular gov UB","Difference LB", "Difference UB") ) )
#(a) Weak and popular government 
s.2a <- sim(m2, x = x.2a, num = c(10000, 100), bootstrap = FALSE)
#Calculate mean
Tab.Mean[1,1] <- mean(as.numeric(s.2a$qi$ev))
#95 percent confidence intervall
est2a.quantiles <- quantile(as.numeric(s.2a$qi$ev), probs= c(.025,.975), na.rm=T) 
Tab.CI[1,1] <- est2a.quantiles[1]
Tab.CI[1,2] <- est2a.quantiles[2]

#(b) Weak and unpopular government 
s.2b <- sim(m2, x = x.2b, num = c(10000, 100), bootstrap = FALSE)
Tab.Mean[1,2] <- mean(as.numeric(s.2b$qi$ev))
est2b.quantiles <- quantile(as.numeric(s.2b$qi$ev), probs= c(.025,.975), na.rm=T)
Tab.CI[1,3] <- est2b.quantiles[1]
Tab.CI[1,4] <- est2b.quantiles[2]

#(ab) Difference between (a) and (b)
s.2ab <- sim(m2, x = x.2a, x1 = x.2b, num = c(10000, 100), bootstrap = FALSE)
Tab.Mean[1,3] <- mean(as.numeric(s.2ab$qi$fd))
est2ab.quantiles <- quantile(as.numeric(s.2ab$qi$fd), probs= c(.025,.975), na.rm=T)
Tab.CI[1,5] <- est2ab.quantiles[1]
Tab.CI[1,6] <- est2ab.quantiles[2]

#(c) Strong and popular government 
s.2c <- sim(m2, x = x.2c, num = c(10000, 100), bootstrap = FALSE)
Tab.Mean[2,1] <- mean(as.numeric(s.2c$qi$ev))
est2c.quantiles <- quantile(as.numeric(s.2c$qi$ev), probs= c(.025,.975), na.rm=T)
Tab.CI[2,1] <- est2c.quantiles[1]
Tab.CI[2,2] <- est2c.quantiles[2]


#(d) Strong and unpopular government 
s.2d <- sim(m2, x = x.2d, num = c(10000, 100), bootstrap = FALSE)  #Predicted prob
Tab.Mean[2,2] <- mean(as.numeric(s.2d$qi$ev))
est2d.quantiles <- quantile(as.numeric(s.2d$qi$ev), probs= c(.025,.975), na.rm=T)
Tab.CI[2,3] <- est2d.quantiles[1]
Tab.CI[2,4] <- est2d.quantiles[2]


#(cd) Difference between (c) and (d)
s.2cd <- sim(m2, x = x.2c, x1 = x.2d, num = c(10000, 100), bootstrap = FALSE)
Tab.Mean[2,3] <- mean(as.numeric(s.2cd$qi$fd))
est2cd.quantiles <- quantile(as.numeric(s.2cd$qi$fd), probs= c(.025,.975), na.rm=T)
Tab.CI[2,5] <- est2cd.quantiles[1]
Tab.CI[2,6] <- est2cd.quantiles[2]


#(ac) Difference between (a) and (c)
s.2ac <- sim(m2, x = x.2a, x1 = x.2c, num = c(10000, 100), bootstrap = FALSE)  
Tab.Mean[3,1] <- mean(as.numeric(s.2ac$qi$fd))
est2ac.quantiles <- quantile(as.numeric(s.2ac$qi$fd), probs= c(.025,.975), na.rm=T)
Tab.CI[3,1] <- est2ac.quantiles[1]
Tab.CI[3,2]  <- est2ac.quantiles[2]


#(bd) Difference between (b) and (d)
s.2bd <- sim(m2, x = x.2b, x1 = x.2d, num = c(10000, 100), bootstrap = FALSE)  
Tab.Mean[3,2] <- mean(as.numeric(s.2bd$qi$fd))
est2bd.quantiles <- quantile(as.numeric(s.2bd$qi$fd), probs= c(.025,.975), na.rm=T)
Tab.CI[3,3] <- est2bd.quantiles[1]
Tab.CI[3,4] <- est2bd.quantiles[2]

#Print results
print(Tab.Mean, digits = 3) #predicted probs
print(Tab.CI, digits = 3) #confidence intervals

#End code Table 2

#############################################################################
#Figure 2: The three-way-interaction plot shows the predicted probability of a dissolution threat varying by the popular support and legislative strength of the government as well as the time remaining until the constitution requires a new election. 
#############################################################################

#Going beyond Table 2, the simulation of predicted probabilities from model 2 in Table 1 now also varies the varible expiration. Otherwise the procedure is as before. 
expirrange <- seq(1,16,1)
#As above
#Public opinion can be low (25th percentile) or high (75th percentile)
po_low <- 0.33
po_hi <- 0.40
#Cabinet seats can be low (25th percentile) or high (75th percentile)
cab_low <- 0.33
cab_hi <- 0.43

#Define layout of plot
par(mfrow=c(1,2)) 

##Start to plot left panel: Government is popular
#Simulation for popular and weak gov
x.2 <- setx(m2, expir = expirrange, po = po_hi, cabseats = cab_low)
s.2 <- sim(m2, x = x.2, num = c(10000, 100), bootstrap = FALSE)
est2.mean <- apply(s.2$qi$ev, MARGIN = 2, FUN=mean)
est2.lb <- apply(s.2$qi$ev,MARGIN=2, quantile, probs= c(.025), na.rm=T)
est2.ub <- apply(s.2$qi$ev,MARGIN=2, quantile, probs= c(.975), na.rm=T)

#Start plot of mean predicted prob
plot(expirrange, est2.mean,
ylim = c(-0.5,1.5),
xlab = "Quarters Until Constitutional End of Term",
ylab = "Pr(Dissolution Threat)",
main = "(a) Government is popular",
pch = 19,
col = 'red',
xlim = c(2,16)
)
abline(h=0)
#Add confidence intervals
segments(expirrange, est2.lb, expirrange, est2.ub, col = 'red')

#Simulation for popular and strong gov
x.2 <- setx(m2, expir = expirrange, po = po_hi, cabseats = cab_hi)
s.2 <- sim(m2, x = x.2, num = c(10000, 100), bootstrap = FALSE)  #Predicted prob
est2.mean <- apply(s.2$qi$ev, MARGIN = 2, FUN=mean)
est2.lb <- apply(s.2$qi$ev,MARGIN=2, quantile, probs= c(.025), na.rm=T)
est2.ub <- apply(s.2$qi$ev,MARGIN=2, quantile, probs= c(.975), na.rm=T)
#Add results to plot
points(expirrange, est2.mean,
pch = 15,
col = 'blue'
)
segments(expirrange, est2.lb, expirrange, est2.ub, col = 'blue')

legend("topright", inset = 0.05, pch = c(19, 15), col = c('red', 'blue'), title = "Seat Share of Government", legend = c("Low", "High"), horiz = TRUE)


##Right panel: Government is unpopular
#Simulation for unpopular and weak gov
x.2 <- setx(m2, expir = expirrange, po = po_low, cabseats = cab_low)
s.2 <- sim(m2, x = x.2, num = c(10000, 100), bootstrap = FALSE)  #Predicted prob
est2.mean <- apply(s.2$qi$ev, MARGIN = 2, FUN=mean)
est2.lb <- apply(s.2$qi$ev,MARGIN=2, quantile, probs= c(.025), na.rm=T)
est2.ub <- apply(s.2$qi$ev,MARGIN=2, quantile, probs= c(.975), na.rm=T)

#Start to plot right panel
plot(expirrange, est2.mean,
ylim = c(-0.5,1.5),
#xlab = "Time (in Quarters) Until Constitutional End of Legislative Term",
xlab = "Quarters Until Constitutional End of Term",
ylab = "Pr(Dissolution Threat) ",
main = "(b) Government is unpopular",
pch = 19,
col = 'red',
xlim = c(2,16)
)
abline(h=0)
segments(expirrange, est2.lb, expirrange, est2.ub, col = 'red')

#Simulation for unpopular and strong gov
x.2 <- setx(m2, expir = expirrange, po = po_low, cabseats = cab_hi)
s.2 <- sim(m2, x = x.2, num = c(10000, 100), bootstrap = FALSE)  #Predicted prob
est2.mean <- apply(s.2$qi$ev, MARGIN = 2, FUN=mean)
est2.lb <- apply(s.2$qi$ev,MARGIN=2, quantile, probs= c(.025), na.rm=T)
est2.ub <- apply(s.2$qi$ev,MARGIN=2, quantile, probs= c(.975), na.rm=T)
points(expirrange, est2.mean,
pch = 15,
col = 'blue'
)
segments(expirrange, est2.lb, expirrange, est2.ub, col = 'blue')

legend("topright", inset = 0.05, pch = c(19, 15), col = c('red', 'blue'), title = "Seat Share of Government", legend = c("Low", "High"), horiz = TRUE)

#End code Figure 2



#############################################################################
#Supporting Information, Table 3: Robustness of the Regression Estimates
#############################################################################

m1 <- zelig(threat ~ po2*cabseats*expir + dpmcp + ssv, model = "logit", data = Data)
m2 <- zelig(threat ~ po2*cabseats*expir + dpmcp + ssv + pmparty + year + threat_lag, model = "logit", data =  Data)
m3 <- zelig(threat ~ po2*cabseats*expir + dpmcp + alternation, model = "logit", data =  Data)
m4 <- zelig(threat ~ po2*cabseats*expir + dpmcp + alternation + pmparty + year + threat_lag, model = "logit", data = Data)
m5 <- zelig(threat ~ po2*cabseats*expir + dcacp, model = "logit", data =  Data)
m6 <- zelig(threat ~ po2*cabseats*expir + dcacp + pmparty + year + threat_lag, model = "logit", data =  Data)
m7 <- zelig(threat ~ po2*cabseats*expir + dpmcp + year + yearsq, model = "logit", data =  Data)
m8 <- zelig(threat ~ po2*cabseats*expir + dpmcp + pmparty + year + yearsq + threat_lag, model = "logit", data =  Data)
m9 <- zelig(threat ~ po2*cabseats*expir + dpmcp, model = "logit", robust = T, data =  Data)
m10 <- zelig(threat ~ po2*cabseats*expir + dpmcp + pmparty + year + threat_lag, model = "logit", , robust = T, data =  Data)
m11 <- zelig(threat ~ po2*cabseats*expir + dpmcp, model = "relogit", bias.correct = T, robust = T, data =  Data)
m12 <- zelig(threat ~ po2*cabseats*expir + dpmcp + pmparty + year + threat_lag, model = "relogit", bias.correct = T, robust = T, data =  Data)

#############################################################################
#Supporting Information, Figure 2
#############################################################################

#Same code as for Figure 2 in paper (above); only replace model input by using Model 11 from abvoe (SI, Table 3)

